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Abstract 

Papillomaviruses (PVs) are widespread pathogens. However, the extent of PV infections in bats remains largely unknown. This work 
represents the first comprehensive study of PVs in Iberian bats. We identified four novel PVs in the mucosa of free-ranging Eptesicus 
serotinus (EserPVI, EserPV2, and EserPV3) and Rhinolophus ferrumequinum (RferPVI) individuals and analyzed their phylogenetic 
relationships within the viral family. We further assessed their prevalence in different populations of £ serotinus and its close relative 
£ isabellinus. Although it is frequent to read that PVs co-evolve with their host, that PVs are highly species-specific, and that PVs do not 
usually recombine, our results suggest otherwise. First, strict virus-host co-evolution is rejected by the existence of five, distantly 
related bat PV lineages and by the lack of congruence between bats and bat PVs phylogenies. Second, the ability of EserPV2 and 
EserPV3 to infect two different bat species (£ serotinus and £ isabellinus) argues against strict host specificity. Finally, the description 
of a second noncoding region in the RferPVI genome reinforces the view of an increased susceptibility to recombination in the E2-L2 
genomic region. These findings prompt the question of whether the prevailing paradigms regarding PVs evolution should be 
reconsidered. 
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Introduction 

Members within the order Chiroptera are extremely successful 
in terms of ecological diversity, accounting for more than one- 
fifth (-1,000) of all extant mammal species (Simmons 2005). 
Bats play a key role in terrestrial ecosystems as pollinators, 
insect controllers, seed dispersers, and reforesters (Kunz 
et al. 201 1). They are also instrumental as vectors of zoonotic 
pathogens, being the reservoir for a number of infectious 
agents capable of crossing species barriers and of infecting 



human and nonhuman hosts (Calisher et al. 2006). Recent 
studies have contributed to enlarge the list of viruses infecting 
bats(Chu etal. 2008; Falcon etal. 2011; Drexler etal. 2012; 
Kurth et al. 2012; Wu et al. 2012), and so far, more than 80 
different viral agents have been identified including the highly 
pathogenic rabies virus (Turner 1975) and related lyssaviruses 
(Samaratunga et al. 1998), Nipah and Hendra viruses 
(Field et al. 2000; Chua et al. 2002), Ebola and Marburg 
viruses (Monath 1999; Leroy et al. 2005), and SARS virus 
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(Li et al. 2005). Moreover, beyond their role as potential res- 
ervoirs of infectious agents, there is growing concern about 
bat disease and mortality. Approximately 25% of the world's 
bat species are threatened with extinction and yet, little is 
known about actual bat pathogens (Rector et al. 2006; 
Muhldorfer et al. 2011). 

Papillomaviridae are a family of small, nonenveloped, 
epitheliotrophic dsDNA viruses. The approximately 8,000 bp 
genome includes an upstream regulatory region, and up to 
eight open reading frames (ORFs) named after their expression 
timing, with an early region encoding for the proteins involved 
in viral replication and immunomodulation and a late region 
encoding for the capsid proteins. Papillomaviruses (PVs) infect 
the skin and mucosa of mammals, but they have also been 
found in birds, turtles, and snakes and probably infect all am- 
niotes (Bravo et al. 2010). Although most PVs cause asymp- 
tomatic infections, some PVs can provoke malignant cell 
transformations. Certain human PVs are responsible for over 
one-third of all infection-associated cancers in humans, includ- 
ing different types of anogenital cancers, head and neck can- 
cers, and skin cancers in genetically susceptible individuals (Zur 
Hausen 2009). Further, neoplastic malignant lesions have also 
been linked to animal PVs in several different hosts: bats 
(RaPV1) (Rector et al. 2006), cats (FcaPV2 and FcaPV3) 
(Lange CE, Tobler K, Markau T, et al. 2009), dogs (CPV1, 
CPV3, and CPV7) (Lange CE, Tobler K, Ackermann M, et al. 
2009), horses (BPV1, BPV2 and EcPV2) (Nasir and Campo 
2008; Scase et al. 2010), rodents (McPV2) (Nafz et al. 
2008), rabbits (SfPV1) (Giri et al. 1985), and sheep (OvPV3) 
(Alberti etal. 2010). 

Reconstructing the evolutionary history of pathogens such 
as PVs requires a taxonomic sampling with a balanced descrip- 
tion of the pathogens' diversity in all possible hosts. However, 
more than half of all known PVs correspond to human PVs 
(PAvE, Papillomavirus Episteme Database, http://pave.niaid. 
nih.gov/, last accessed December, 2013). Considering that 
there are over 23,000 amniote species serving as potential 
hosts (Bravo et al. 2010; International Union for 
Conservation of Nature 2013, http://vvww.iucnredlist.org/, 
last accessed December, 2013), we have little insight of 
nonhuman PV diversity. The development of a comprehensive 
and unbiased scenario for PVs evolution demands the consid- 
eration of poorly studied groups such as bats. To date, only 
five different bat PVs have been fully sequenced. MschPVI, 
MschPV2, and MrPV1 were isolated, respectively, from oro- 
pharyngeal and/or anal swabs from healthy free-ranging 
Miniopterus schreibersii and Myotis ricketti individuals (Tse 
et al. 2012; Wu et al. 2012), while EhelPVI was retrieved 
from hair bulbs from a healthy captive Eidolon helvum indi- 
vidual (Garcia-Perez et al. 2013). RaPV1 was recovered from a 
basosquamous carcinoma on the wing of a Rousettus aegyp- 
tiacus individual (Rector et al. 2006). Partial sequences of the 
E1 and L1 genes have also been retrieved from hair bulbs from 
a healthy captive Pteropus giganteus (Garcia-Perez et al. 



201 3). Further, Baker et al. (201 3) conducted a metagenomic 
study on African Ei. helvum bats, resulting in the identification 
of hundreds of short PV-like sequences, isolated from throat, 
lung, and urine samples. 

In this study, we have systematically surveyed the presence 
of PVs in oropharyngeal swabs from 22 out of the 31 extant 
Iberian bat species (VV.AA. 2010). We have identified and 
completely sequenced the complete genomes of four novel 
PVs isolated from two different bat species, namely Eptesicus 
serotinus and Rhinolophus ferrumequinum. Moreover, we 
have assessed the prevalence of the novel PVs in Iberian col- 
onies of E. serotinus and the morphologically very similar, 
closely related species E. isabellinus (Juste et al. 2013). 

Material and Methods 

Ethics Statement 

All persons responsible for samples collection were qualified 
and experienced bat researchers who had bat capture and 
sampling permits issued by the competent environmental au- 
thority of their study regions as follows: Direccion General de 
Gestion del Medio Natural, Consejeria de Medio Ambiente, 
Junta de Andalucia, Spain (code #201230E040); Direccion 
General de Montes y Espacios Naturales, Consejeria de 
Agricultura, Junta de Comunidades de Castilla-La Mancha, 
Spain (code #DGMEN/SEN_avp_12083_aut); Direccion 
General del Medio Natural, Consejeria de Fomento y Medio 
Ambiente, Junta de Castilla y Leon, Spain (code #EP/CYL/201/ 
2012); Direccion General de Medio Ambiente, Consejeria de 
Agricultura, Desarrollo Rural y Medio Ambiente y Energia, 
Gobierno de Extremadura, Spain (code #CN009/12/ACA). 
The sampling protocol was approved by the Bioethical 
and Animal Welfare Committee (CEBA-EBD) of the Estacion 
Biologica de Doriana (EBD-CSIC), study code #CEBA- 
EBD_11_30, adhering to the guidelines in the RD1201/2005 
on the protection of animals used for experimentation and 
other scientific purposes. 

Sample Collection and Nucleic Acid Extraction 

A total of 44 oropharyngeal swabs were taken from free- 
ranging bats belonging to 22 different species through 
Spain during 2002-2008. An exhaustive list of the samples 
is provided in supplementary table S1, Supplementary 
Material online. Bats identification was mainly based on mor- 
phological characters. For sibling or morphologically cryptic 
complexes, the identification was based on the sequencing 
of a diagnostic mitochondrial fragment (Cytb -500 bp) follow- 
ing the protocol established by Ibariez et al. (2006). Sample 
collection and DNA extraction were performed as described 
previously (Echevarria et al. 2001). Several contention mea- 
surements were undertaken in order to prevent contamina- 
tion: all pipetting steps were performed under safety hood 
cabinets; DNA extraction and subsequent DNA amplification 
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were performed in different locations (i.e., Madrid and 
Barcelona); all pre- and post-PCR manipulations were 
performed in different facilities. 

Viral Genome Amplification, Sequencing, and Cloning 

Rolling circle amplification (RCA) was performed to generate 
full-length PV genomes following an optimized protocol 
(Schulz et al. 2009). PV presence was tested using the PV 
degenerated FAP (Forslund et al. 1999) and CP primers 
(Iftner et al. 2003). The FAP and CP primers amplify DNA frag- 
ments of approximately 450 bp of the L1 and E1 genes, re- 
spectively. Complete viral genome amplification, sequencing, 
and cloning were performed as described previously (Garcia- 
Perez et al. 201 3). The complete genomes were sequenced by 
primer walking using forward and reverse primers, and each 
position was read at least twice in each direction. Overlapping 
fragments were generated and used to clone each of the four 
novel PVs using the CloneJet PCR Cloning Kit (Fermentas) and 
DH5a cells (Invitrogen). Two clones per overlapping fragment 
were selected and sequenced to further confirm the PVs 
genome sequence. Novel PVs were designated EserPVI, 
EserPV2, EserPV3, and RferPVI, following suggestions to 
avoid naming ambiguities (Garcia-Perez et al. 2013). 

Genomic and Protein Sequence Annotation 

The ORFs encoded in EserPVI , EserPV2, EserPV3, and RferPVI 
were identified with the ORF Finder tool on the NCBI server. 
They were confirmed by comparison with a nonredundant 
protein sequences database in GenBank through the 
BLASTP server. The MEME algorithm (Bailey et al. 2009) was 
used to identify putative E2 binding site (E2BS) sequence pat- 
terns occurring in Lambda + MuPVs and previously described 
bat PVs. Pairwise identities and similarity values for nucleotide 
and protein sequences were calculated using the EMBOSS 
Needle software. 

Phylogenetic Analysis 

The data set used to examine the phylogenetic relationships of 
the four novel viruses comprised 143 PVs that covered their 
currently known diversity in terms of hosts. A comprehensive 
description of the PVs employed in this study is provided in 
supplementary table S2, Supplementary Material online. For 
this selection, the E1, E2, L2, and L1 genes were analyzed. 

Amino acid alignments were constructed individually 
for each protein (E1, E2, L2, and L1) with MUSCLE (Edgar 
2004), filtered with GBLOCKS (Castresana 2000), and conca- 
tenated. Maximum likelihood phylogenetic inference for the 
E1-E2-L2-L1 combination was conducted with RAxMLv7.2.8 
(Stamatakis 2006), using the GTR + T4 substitution model for 
nucleotide alignments and the LG substitution model 
for amino acid alignments (Gottschling, Goker et al. 2011). 
Phylogenetic analyses were calculated considering the 
corresponding partitions (12 for nucleotide alignments, 



corresponding to one per codon position per gene, and 4 
for amino acid alignments), with individual per partition 
branch length optimization and running 1 ,000 bootstrap rep- 
licates. Additional phylogenetic analyses were performed as 
described above separately for the E1-E2 and L2-L1 combi- 
nations to investigate possible topological incongruences. 

Individual phylogenies for the E6, E7, E1, E2, L2, and L1 
genes were constructed exclusively for the Lambda + MuPVs 
and used to create a supernetwork using SPLITSTREE v4 
(Huson and Bryant 2006). Bayesian inference was performed 
on the same alignments with PHYLOBAYES v3.3 (Lartillot et al. 
2009) applying the GTR + T4 substitution model for nucleo- 
tide alignments and the LG substitution model for amino acid 
alignments, removing constant sites, running two indepen- 
dent chains, and checking for convergence comparing dis- 
crepancies among partitions. The convergence criterion 
between chains was that the largest discrepancy observed 
across all bipartitions should be below 0.1. Trees were 
rooted a posteriori using the sequences of PVs found in 
birds and in turtles. 

A phylogeny including the bat species from which com- 
plete PV genome sequences have been retrieved was con- 
structed using a concatenation of partial sequences from 
three mitochondrial markers, namely Cytochrome Oxidase I 
(COI, 592 bp), Cytochrome b (Cytb, 754 bp), and NADH de- 
hydrogenase (ND1, 51 8 bp), and the nuclear gene recombi- 
nation activating gene (RAG2, 759 bp). Bayesian inference 
was performed with BEAST applying the GTR substitution 
model, using four partitions (one per gene), running two in- 
dependent chains and empirically estimating the gamma dis- 
tribution parameter. The first 1 00,000 trees were discarded as 
burn-in, and posterior probabilities were calculated by sam- 
pling 3 x 10 6 generations every 300 trees. 

All viral contig sequences obtained by Baker et al. (2013) 
were screened for PV origin using TBIastX searches, further 
identifying for each contig the coding frames and sequences 
and eliminating premature stop codons. The final curated 
amino acid sequences are provided in supplementary 
table S3, Supplementary Material online. An evolutionary 
placement algorithm (EPA) (Berger and Stamatakis 2011) 
was applied to introduce the sequences into the previously 
constructed well-resolved phylogeny as described (Mengual- 
Chuliaetal. 2012). A total of 381 sequences (E1 n=151,E2 
n = 40, L2 n = 55, and L1 n= 135), together with the partial 
E1 and L1 sequences of PgigPVI, were assigned a phyloge- 
netic position. The same strategy was used to infer the phy- 
logenetic positions of the partial E1 (n = 5) and L1 (n = 9) 
sequences obtained in this study. 

Prevalence of EserPVI, EserPV2, and EserPV3 in Iberian 
Bat Populations 

Prevalence was investigated by PCR screening of 267 addi- 
tional samples, including oropharyngeal swabs (n = 78), 
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anogenital swabs (n = 85), and hair bulbs (n= 104) obtained 
from 1 06 £ serotinus (n = 33) and £ isabellinus (n = 73) indi- 
viduals. Samples were collected from central Spain during 
2012. A detailed summary of the samples is provided in sup- 
plementary table S4, Supplementary Material online. DNA ex- 
traction from swabs was performed as described above. DNA 
extraction from hair bulbs was performed following a protocol 
for the isolation of genomic DNA from tissues (Qiagen). Full- 
length PV genomes were amplified using RCA, and RCA prod- 
ucts were used as template for the PCR reactions. Primer se- 
quences, gene targets, amplicon sizes, and PCR conditions are 
specified in supplementary table S5, Supplementary Material 
online. Presence of other unknown PVs was further studied 
using the PV-specific FAP and the CP primers, as described 
(Mengual-Chulia et al. 2012). Obtained amplicons were also 
sequenced in both strands using the same amplification pri- 
mers and analyzed through a BlastX search to confirm their PV 
origin. 

Results 

Genomic Organization and Sequence Similarity to 
Other PVs 

We have characterized the complete genomes of four novel 
PVs isolated from mucosal swabs of two free-ranging 
Iberian £ serotinus (EserPVI , EserPVZ, and EserPV3) and one 
R. ferrumequinum (RferPVI) individuals. 

The genomes of EserPVI, EserPV2, EserPV3, and RferPVI 
consisted of 7,668, 7,574, 7,71 1, and 8,249 bp, respectively, 
and are available in GenBank under the following accession 
numbers: KC858263, KC858264, KC858265, and 
KC858266. The four genomes presented the typical PV 
ORFs coding for five early proteins (E6, E7, E1, E2, and E4) 
and two late proteins (L2 and L1) located on the same coding 
strand. Putative E4 ORFs were identified nested within E2, as 
they contained rich proline stretches that characterize E4 pro- 
teins. The four PVs presented a noncoding region (NCR1) 
spanning between the stop codon of L1 and the start 
codon of E6. Additionally, RferPVI contained a second 
noncoding region between the early and late regions. The 
genomic arrangement of these four PVs and the location of 
some common amino acid motifs and regulatory elements are 
depicted in supplementary figure S1, Supplementary Material 
online. A detailed description of the genomes, including 
the precise location of ORFs, amino acid motifs, and regula- 
tory elements, is provided in supplementary table S6, 
Supplementary Material online. 

Two classical zinc-binding domains (CX 2 CX 2 9CX 2 C) sepa- 
rated by 36 amino acids were located in all E6 oncoproteins. 
The same motif was found in the E7 oncoproteins of RferPVI 
and EserPV2 and slightly modified (CX2CX30CX2C) in the E7 
proteins of EserPVI and EserPV3 and in those of EhelPVI, 
MschPV2, and RaPV1 . The E7 oncoproteins of EserPV2 and 



RferPVI contained a potential pRB-binding motif (LXCXE), 
absent in the E7 proteins of EserPVI, EserPV3, and MschPV2. 
The ATP-binding domain (GX 4 GKS) appeared in all E1 proteins. 
A leucine zipper domain (LX 6 LX 6 LX 6 L) was found in the E2 
proteins of EserPVI and EserPV3 but was absent in EserPV2, 
RferPVI EhelPVI, MschPVI, and MschPV2 (supplementary 
table S6, Supplementary Material online). 

Polyadenylation signals (AATAAA) for the early and late 
transcripts were found at the beginning of the L2 gene and 
in the NCR1, respectively. TATA boxes (TATAAA) of the E6 
promoter were located close to the 3' end of the NCR1. 
Canonical (ACCG-N 4 -CGGT) and noncanonical putative E2- 
binding sites (E2BS) were detected in the NCR1 and also in the 
L2 genes of the novel bat PVs. The canonical and noncanon- 
ical E2BS patterns occurring in the NCR1 and L2 genes of 
Lambda + MuPVs, and the E2BS pattern occurring in the 
NCR1 of bat PVs, are illustrated in supplementary figures S2 
and S3, Supplementary Material online. The presence of pu- 
tative E2BS had not been previously reported within the L2 
genes of PVs, and their biological significance is unknown. 

Sequence similarities among all bat PVs described to date 
were investigated by pairwise alignments of the E6, E7, E1 , E2, 
L2, and L1 genes and of their respective proteins (supplemen- 
tary table S7, Supplementary Material online). Sequence sim- 
ilarities for the L1 genes and proteins were also studied within 
the Lambda + MuPVs crown group (supplementary table S8, 
Supplementary Material online). EserPVI and EserPV3 showed 
the highest similarity to each other while RferPVI shared its 
highest similarity with both EserPVI and EserPV3. EserPVZ was 
most similar to MschPVI . Low similarity percentages were 
found between the novel PVs and the bat PVs MrPV1 and 
RaPV1 . Considering the ICTV guidelines for delineating PV 
taxonomy, based exclusively on nucleotide identity on the L1 
gene (de Villiers et al. 2004), the novel PVs here described 
RferPVI, EserPVI, and EserPV3 could belong together into a 
single novel genus within the Lambda + MuPVs crown group, 
while EserPV2 and MschPVI belong together into a different 
genus, branching close to the origin of the four PV crown 
groups (supplementary tables S7 and S8, Supplementary 
Material online). 

Phylogenetic Analysis 

Phylogenetic trees were calculated using maximum likelihood 
and Bayesian approaches, both at the nucleotide and at the 
amino acid levels. Details of the different alignments and par- 
titions considered are provided in supplementary table S9, 
Supplementary Material online. Topogical incongruences be- 
tween early (E1-E2) and late (L2-L1) phylogenies were not 
observed for the novel bat PVs. The E1-E2-L2-L1 concatena- 
tion was therefore used for subsequent analyses. All recon- 
structed phylogenetic trees are available in the supplementary 
material, Supplementary Material online, as suggested (Drew 
et al. 2013). The E1-E2-L2-L1 gene combination generated 
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well-supported phylogenies (fig. 1) in which PVs segregated 
into four major groups, as previously described: 
Alpha + OmikronPVs (infecting Artiodactyla, Carnivora, 
Cetacea, Chiroptera, and Primates), Beta + XiPVs (infecting 
Artiodactyla, Carnivora, Eulipotyphla, Primates, and 
Rodentia), Delta + ZetaPVs (infecting Artiodactyla and 
Perissodactyla), and Lamba + MuPVs (infecting Carnivora, 
Chiroptera, Lagomorpha, Primates, and Rodentia). Bat PVs 
did not constitute a monophyletic group. Instead, phyloge- 
netic analyses revealed the existence of at least five different 
bat PV lineages: 1) RaPV1 showed an uncertain position basal 
to the tree; 2) EserPV2 and MschPVI appeared as sister taxa, 
close to the basal branching events of the four crown 
groups; 3) MrPV1 was the sister taxon of the UmPV1, confi- 
dently nested within the Alpha + OmikronPVs crown group; 4) 
EserPVI, EserPV3, and RferPVI constituted a monophyletic 
clade within the Lambda + MuPVs crown group; 5) EhelPVI 
and MschPV2 showed an uncertain phylogenetic position 
within the Lambda + MuPVs crown group (fig. 1 and supple- 
mentary fig. S4, Supplementary Material online). The phylo- 
genetic network reconstructed for Lambda + MuPVs 
confirmed the closer evolutionary relationships among 
EserPVI, EserPV3, and RferPVI but could not resolve the am- 
biguity regarding EhelPVI and MschPV2 (fig. 2). Bat PV-like 
sequences obtained by Baker et al. (2013) in their metage- 
nomic study of the Ei. helvum virome appeared scattered 
throughout the PV phylogeny and belonged to all different 
crown groups (supplementary fig. S5, Supplementary Material 
online). 

Prevalence of EserPVI, EserPV2, and EserPV3 in Iberian 
Bat Populations 

Anogenital, oropharyngeal, and hair bulb samples from 33 
E. serotinus and 73 E. isabellinus individuals were taken 
from seven Iberian bat colonies. All samples were tested for 
the presence of EserPVI, EserPV2, and EserPV3 DNA using 
specific primers. The amplicons and regions targeted for 
each virus are detailed in supplementary table S5, 
Supplementary Material online. In total, 10% (8/78) of oro- 
pharyngeal samples, 6% (5/85) of anogenital samples, and 
none of the hair bulbs tested positive for DNA of any of 
these viruses. Our results highlight the essentially mucosal tro- 
pism of these three PVs, also initially retrieved from oropha- 
ryngeal swabs. Regarding species specificity, 30% (10/33) 
among the E. serotinus individuals and 5% (3/62) among 
the E. isabellinus individuals tested positive for DNA of any 
of these three viruses (fig. 3 and supplementary table S10, 
Supplementary Material online). All amplified sequences cor- 
responded exactly to the one isolated originally, independent 
of the host species. Additionally, presence of other unknown 
PVs on the same samples was studied using the PV-specific 
FAP and CP primers. In total, 21 % (16/78) of oropharyngeal 
samples and 0% of anogenital samples tested positive for the 



presence of DNA of other PVs (fig. 3 and supplementary 
table S10, Supplementary Material online). The accession 
numbers for the partial L1 and E1 sequences amplified, re- 
spectively, by the FAP and CP primers are provided in supple- 
mentary table S11, Supplementary Material online, and 
pairwise nucleotide identities between these partial sequences 
and their closest relatives are provided in supplementary 
table S12, Supplementary Material online. The alignment of 
the novel fragments amplified with the FAP and CP primers is 
provided in supplementary table S1 3, Supplementary Material 
online. Phylogenetic analysis of these short PV sequences re- 
vealed that they were most closely related to the novel bat PVs 
reported here and to the Sigma-PV EdPV1 and the Nu-PV 
HPV41 (supplementary fig. S6, Supplementary Material 
online). 

Discussion 

There is an intrinsic value in the explicit introduction of the 
study of pathogen evolutionary history into molecular medi- 
cine research. Understanding the phylogenetic relationships 
and mechanisms driving PVs diversification will provide 
better insight into viral evolution and virus-host interactions. 
Three assumptions have traditionally dominated PV research: 
that PVs have co-evolved with their hosts (Van Ranst et al. 
1995); that PVs are highly species-specific (Van Ranst et al. 
1995; Halpern 2000; Bernard et al. 2006); and that there are 
minimal virus-virus interactions and therefore events such as 
recombination are rare (Plummer et al. 2011). Although all 
three dogmas have been recurrently challenged by recent re- 
sults (Gottschling et al. 2007; Woolford et al. 2007; 
Gottschling, Goker et al. 2011; Sakakibara et al. 2013), it is 
still a commonplace to assume their validity and to address 
research on human PVs as if they were a monophyletic, dis- 
tinct entity from animal PVs (Brody 2012; Bernard 2013; 
de Villiers 2013). 

Virus-Host Co-evolution Is Not a Major Determinant of 
Mammalian PV Evolution 

The PV-host co-evolution hypothesis has relied on the gross 
phylogenetic correspondence for certain PVs and for their 
mammalian hosts: Alpha- and BetaPVs infecting primates, 
DeltaPVs infecting ruminants, LambdaPVs infecting carni- 
vores, and PiPVs infecting rodents. However, a systematic to- 
pology analysis reveals the absence of congruence between 
these PVs and their host phylogenies, without a single exam- 
ple of identical tree topologies for both PVs and hosts (Chan 
et al. 1997; Antonsson and Hansson 2002; Gottschling, Bravo 
et al. 201 1; Gottschling, Goker et al. 201 1). The contribution 
of the different mechanisms to the evolution of PVs can be 
quantified, and virus-host co-evolution explains only around 
30% of all events needed to invoke reconciliation of PVs and 
hosts phylogenetic trees (Gottschling, Goker et al. 201 1). This 
value should nevertheless be taken with caution, given the 
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Fig. 1. — Bayesian amino acid phylogenetic reconstruction for the E1-E2-L2-L1 concatenation. Branch lengths are drawn to scale, with the scale 
bar indicating the evolutionary distance in substitutions per site. Numbers above the branches indicate Bayesian posterior probabilities and ML bootstrap 
support values. Maximum support values are indicated with an asterisk (*), while values below 0.50 and 50 are indicated with a dash (-). Color code 
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limited taxon of animal PVs and the focus on human PVs. The 
lack of global congruence between PVs and host trees is also 
reflected on the existence of several polyphyletic lineages for 
PVs infecting the same host. This is true for PVs isolated from 
humans, chimpanzees, gorillas, macaques, rodents, dogs, 
cats, cattle, sheep, horses, dolphins, and porpoises (fig. 1 
and supplementary fig. S4, Supplementary Material online), 
which belong to different PV genera and appear scattered 
throughout the PV phylogeny in a highly polyphyletic pattern 
(Bravo et al. 2010; Gottschling, Goker et al. 201 1). 

The description of the novel bat PVs communicated here 
provides further evidence to reject the hypothesis of exclu- 
sively PV-host co-evolution. Bats monophyly is widely ac- 
cepted and supported both by morphological (Simmons 
1994) and molecular studies (Murphy et al. 2001). The exis- 
tence of several polyphyletic bat PVs lineages rejects the 
hypothesis of co-evolution between PVs and bats at a global 
scale. Likewise, co-evolution as a more recent event seems 
also unlikely. Under a strictly co-evolution pattern, one 
would expect pteropid and rhinolophid PVs (RaPV1, 
EhelPVI, and RferPVI) to be monophyletic, while PVs infect- 
ing the remaining bats should form a different clade. One 
would also expect serotine and myotid PVs to be more closely 
related and both to be sister taxa to the miniopterid PVs 
(fig. 4). Instead, PVs infecting pteropids, rhilonophids, and 
vespertilionids are found intermingled in several, only distantly 
related, PV lineages (fig. 1). In addition to the four novel 
genomes here communicated, we have identified fragments 
of ten possible novel PVs from E. serotinus and four from 
E. isabellinus (supplementary tables S4 and S11, 
Supplementary Material online). Furthermore, the phyloge- 
netic positions assigned here to other short PV-like sequences 
isolated from bats (Baker et al. 201 3) suggest that the number 
of different bat PV lineages might be even larger (supplemen- 
tary fig. S5, Supplementary Material online). The performed 
analysis placed the metagenomic sequences in close proximity 
to the novel bat PVs described, but strikingly, some taxa 
appeared nested within all PVs crown groups. The appearance 
at the same position of sequences retrieved from different 
PV genes supports these different locations. Non-assembled 
sequences must however be regarded preliminary, and 
cloned, complete genome sequences would be required 
to fully confirm the relationships of the putative novel 
viruses. Yet, it seems likely that bats can be the hosts to a 
plethora of different PVs belonging to all different crown 
groups. 



Remarkably, all isolates containing EserPV2 DNA corre- 
sponded to coinfections between this virus and a second PV 
(supplementary fig. S6 and table S4, Supplementary Material 
online). Interestingly, serotine bats from one of the sampled 
colonies (Casatejada) showed abundant ulcer-like wounds in 
their wings (data not shown). This was also the colony with 
the highest prevalence of PV infection, with 50% of the indi- 
viduals testing positive for any PV (supplementary tables S4 
and S10, Supplementary Material online). Further research is 
being carried out to establish the etiology of these lesions, 
unlikely to be associated to PV infections. The PV-specific 
FAP and CP primers failed to detect EserPVI, EserPV2, and 
EserPV3 in the screened samples. This fact suggests, in agree- 
ment with previous studies, that the efficiency of the universal 
consensus FAP and CP primers is limited and that they miss a 
considerable number of infections. It is therefore very likely 
that the prevalence of PV infection is larger than estimated 
here. 

No Strict Species Specificity in PVs 

PVs are usually considered as highly species-specific (Bernard 
et al. 2010; Bernard 2013). PVs are named after the host in 
which they were first isolated, following the ICTV code 
(Bernard et al. 2010). All subsequent analyses considering 
host-virus specificity commonly follow this naming conven- 
tion as a hypothesis to assume host specificity or at least 
host preference. However, the results communicated here 
show that similar viruses can infect different hosts, which 
sums up to the several examples of heterologous PV infections 
between distantly related hosts (Munday et al. 2007; van Dyk 
et al. 2009; Munday and Knight 2010; Gottschling, Bravo 
et al. 2011). We initially isolated EserPVI, EserPV2, and 
EserPV3 from oropharyngeal swabs from two different 
E. serotinus individuals. Notably, the results of the screening 
performed in bat samples collected from seven different 
colonies revealed the presence of DNA from EserPV2 and 
EserPV3 not only in E. serotinus but also in E. isabellinus indi- 
viduals. These are two morphologically cryptic but molecularly 
highly differentiated species (above 16% divergence in the 
cytochrome b gene) within the serotinus species group with 
different geographic distributions. E. serotinus is widely distrib- 
uted throughout Europe, while E. isabellinus shows a patchy 
distribution restricted to the southern half of the Iberian 
Peninsula and North Africa (Ibahez et al. 2006; Juste et al. 
2013). The growing number of PV cross-infections descrip- 
tions suggests that certain PVs might exhibit a broader host 



Fig. 1. — Continued 

highlights the four PV crown groups: red, Alpha + OmikronPVs; green, Beta + XiPVs; blue, Delta + ZetaPVs; ochre, Lambda + MuPVs. Viruses whose detailed 
phylogenetic relationships could not be disentangled are labeled in black. Silhouettes represent the infected hosts. Taxonomic classification of both hosts 
(host order) and viruses (PV genera) are included. Gray dots highlight the five lineages encompassing bat PVs. Branches corresponding to clades or PVs that 
contain an E2-L2 region and may thus reflect individual recombination events are highlighted with a black star. The novel bat PVs described here are 
highlighted with black arrows. 
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Fig. 2. — Lambda + MuPVs supernetwork. The network was constructed using the best-known maximum-likelihood trees of each individual nucleotide 
PV gene (E6, E7, E1, E2, L2, and LI). Color code represents the different orders of the hosts. Specific PVs tropisms and outcome of the corresponding 
infections are indicated in the inset. PV genera are specified in gray. 
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Fig. 3— Prevalence of EserPVI , EserPV2, EserPV3, and other PVs DNA 
in the screened samples recovered from seven different Iberian £ serotinus 
and F. isabellinus colonies. 
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Fig. 4. — Tanglegram linking the phylograms of bat PVs and their 
hosts. Congruence between the phylogenetic relationships among host 
bat species after Bayesian inference (left) and bat PVs after maximum 
likelihood inference (right). Please note that the PV taxa depicted here 
are representatives of highly polyphyletic PV crown groups. Color code 
for the lines linking both phylograms corresponds to colors used in figure 2 
for the different PV crown groups. 
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range and/or that host switches between different species 
might occur more frequently than initially suspected. This 
potential to infect different hosts could help explain the 
highly polyphyletic and/or paraphyletic pattern that many 
PVs display. 

Presence of a Second NCR in RferPVI 

During evolution, and independently for different PV clades, 
the E2-L2 intergenic region has often behaved as a target for 
recombination: during viral integration in cancers (Doorbar 
et al. 2012); early in the evolution of AlphaPVs (Bravo and 
Alonso 2004; Varsani et al. 2006; Angulo and Carvajal- 
Rodriguez 2007); between PVs belonging to different 
genera (Gottschling, Bravo et al. 2011; Gottschling, Goker 
et al. 2011; Robles-Sikisaka et al. 2012); and even between 
members of two viral families, Papillomaviridae and 
Polyomaviridae (Woolford et al. 2007). The E2-L2 region 
can also accommodate both coding and noncoding genomic 
segments, which may have gained access to the PV genomes 
through recombination events with hitherto nonidentified 
donors. On the one hand, Alpha- and DeltaPVs encode in 
their E2-L2 region for different nonevolutionarily related E5 
proteins (Bravo and Alonso 2004). On the other hand, 
LambdaPVs, EePV1, isolated from an European hedgehog, 
and the RferPVI reported here, all incorporate a large noncod- 
ing region (NCR2) in the intergenic E2-L2 region. No regula- 
tory, promoter, or coding elements can be identified in these 
NCR2, and the presence and conservation of such a long seg- 
ment of around 1 kb, i.e., 15% of the viral genome, is puz- 
zling. All these putative recombination events in the intergenic 
E2-L2 region of the PVs have likely occurred as individual, 
independent events. The alternative hypothesis of all clades 
containing elements between E2 and L2 to be monophyletic is 
rejected (Shimodaira-Hasegawa test, P value <0.01) when 
compared with nonconstrained trees. Our interpretation 
implies thus that at least five independent recombination 
events (identified in fig. 1) may have occurred in between 
the E2 and L2 genes throughout the evolutionary history 
of PVs. 

Conclusion 

The description here reported of four novel PVs and the as- 
sessment of their prevalence in Iberian bat populations rein- 
force the notion that the evolutionary dynamics of PVs are 
complex. Our results strongly point toward multiple forces 
as drivers of PVs evolution, including co-evolution, adaptive 
radiation, broad host range, host switch, and recombination. 
However, as long as our knowledge of PVs, diversity is biased 
toward a specific host (humans) and thus unbalanced, sound 
conclusions about PVs evolution cannot be reached. The 
number of nonhuman PVs has increased in recent years (cur- 
rently 114 known nonhuman PVs isolated from 55 different 
species). Yet, taxonomic sampling is insufficient and clades 



such as Afrotheria or Xenarthra remain largely unexplored. 
A systematic sampling would eventually allow the develop- 
ment of a comprehensive evolutionary framework reconciling 
the biology, epidemiology, and genomic structure of PVs. Only 
with this scenario we will be able to provide an answer to 
many of the still unsolved questions such as the presence of 
signatures in different PVs that account for their oncogenic 
potential, their cell tropism, or their host range. 

Supplementary Material 

Supplementary figures S1-S6 and tables S1-S13 are available 
at Genome Biology and Evolution online (http://www.gbe. 
oxfordjournals.org/). 
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